function [N_prefs, N_alphas, alpha_list, start_H] = make_Q_c2( param )
% Generates the objective matrix Q given payoff parameters.



%% SETUP

% Global variables:

global samp
global Shocks
global N_types l_t first_wproj can_add Add_cell
global N N_Z feasible isolated residual1 residual2 Z_names
global Y_t nco_t dif_t 

% Preference parameters:

al = param(1);  % expected impact
c0 = param(2);  % constant cost
cN = param(3);  % additional coauthors
cD = param(4);  % unfamiliar topic

% Increasing cost of links:
% (calibrated to offset mean of sorted shocks)

c1 = -0.6;

% Min/max values for combined categories of types:

aY_max = max(max(al*Y_t) );

c_min = c0 + min( [0, cN*1, cN*2, cN*3] ) + min(0, cD);



%% PREFERENCE CLASSES

% Initialize variables:

N_prefs = 0;

clear global N_alphas
global N_alphas

N_alphas = 0;

% Loop over individual characteristics:

tic
for a = 1:size(Z_names,1)
	z = Z_names(a,:);
	N_z = N_Z.(z);
    % Lists of types:
	feas0_z = int32( feasible.(z) );
	isola_z = int32( isolated.(z) );  % (not used)
	resid1_z = int32( residual1.(z) );
	resid2_z = int32( residual2.(z) );
    % Actual types with links:
	feas1_z = feas0_z(feas0_z > max(resid2_z) );
    % Actual or residual types with links:
    feas2_z = [resid1_z ; resid2_z ; feas1_z];
	% Common utility for actual types:
	ubaseA = al*Y_t(feas1_z,:) - ...
		c0 - cN*nco_t(feas1_z,:) - cD*dif_t(feas1_z,:);
	% Marginal project for actual types:
	umargA = min(ubaseA,[],2);
	% Marginal project for residual categories (max possible):
	umargR1 = double(resid1_z*0)+ aY_max - c_min;
	% Common utility for combined categories (uses max expected impact):
	ubaseR2 = aY_max - ...
		c0 - cN*nco_t(resid2_z,:) - cD*dif_t(resid2_z,:);
	% Marginal project for specific residual types:
	umargR2 = min(ubaseR2,[],2);
    % Marginal project for all sets of types:
    umargB = [umargR1 ; umargR2 ; umargA];
	% Subtract cost shock (and increasing cost of links), is net positive:
	ok_i = repmat(umargB',length(Shocks),1) - repmat(c1*l_t(feas2_z)',length(Shocks),1) - Shocks(:,l_t(feas2_z)') > 0;
	% Unique preference classes (across draws of shocks):
	[ok_z,~,ids] = unique(ok_i,'rows');
    % Include isolated type:
	ok_z = [true(size(ok_z,1), 1) ok_z];
	% Compute probabilities:
	temp = tabulate(ids);
	muZpH_z = temp(:,3) / 100 * N_z/N;
	% Store results:
	N_prefs = N_prefs + size(ok_z, 1);
	[cl,rw] = find(ok_z');
	N_alphas = N_alphas + size(cl, 1);
	assignin('caller', ['H_' z], [rw feas0_z(cl)]);
	assignin('caller', ['muZpH_' z], muZpH_z);
end
disp('Preference classes:')
toc



%% ALLOCATION PARAMETERS

tic

% Initialize variables:

alpha_list = zeros(1,N_alphas,'int32');

next = 0;

% Make combined list of allocation parameters:
% (vector of types contained in each pref. class)

for a = 1:size(Z_names,1)
	z = Z_names(a,:);
	H_Z = evalin('caller', ['H_' z]);
	len = size(H_Z,1);
	alpha_list( next + (1:len) ) = H_Z(:,2);
	next = next + len;
end

% First element of each preference class in above list:
% (an isolated type)

start_H = find(alpha_list < first_wproj);

% Positions of each network type in list of allocation parameters;

tofa_idx = cell(N_types,1);

for t = find(can_add)'
	tofa_idx{t} = find(alpha_list == t);
end

disp('Allocation parameters:')

toc



%% S MATRIX and muZpH_vec
% (S [N_alphas x N_types]: indicates desired alter types for each alpha)

% Approximate number of nonzero elements:

N_can_add = sum( can_add( alpha_list) );
N_nz = round( N_can_add * N_types / 10 );
%%%NOTE: Dividing by 10 b/c not all types add to all other types.

% Initialize variables:

S_cells = zeros(2,N_nz,'uint32');

clear global muZpH_vec
global muZpH_vec

muZpH_vec = zeros(N_alphas,1);

% Initialize counters:

next = 0;
at_idx = 0;
at0_idx = 0;

% Loop to find desired alter types for each type in each preference class
% and make vector of preference class probabilities:

tic
for a = 1:size(Z_names,1)
	z = Z_names(a,:);
    % All pref. classes for this Z:
	H_Z = evalin('caller', ['H_' z]);
	muZpH_z = evalin('caller', ['muZpH_' z]);
    % Loop over these pref. classes:
	for h = 1:max(H_Z(:,1) )
        % One pref. class:
		H = H_Z(H_Z(:,1) == h, 2);
        % Repeat pref. class probability by its number of elements:
		muZpH_vec( at0_idx + (1:size(H,1)) ) = repmat(muZpH_z(h),1, size(H,1) );
        can_add_idx = find(can_add(H));
        % Loop over types that can add projects:
        for b = can_add_idx'
			at_idx = at0_idx + b;
            t = H(b);
            % Extract desired alter types:
            temp = Add_cell(t,H);
			s_want = unique( cell2mat( temp(~cellfun('isempty', temp) )' ) );
            % Add them to list of index positions for S matrix:
            if ~isempty(s_want)
				len = length(s_want);
                S_cells(1, next + (1:length(s_want)) ) = repmat(at_idx,1,len);
				S_cells(2, next + (1:len) ) = s_want;
				next = next + len;
            end
        end	
        at0_idx = at0_idx + size(H,1);
	end
end
disp('S matrix cells (and muZpH_vec):')
toc

% Make sparse matrix from list of index positions:

tic
S = sparse(double(S_cells(1,1:next)),double(S_cells(2,1:next)), true, N_alphas,N_types);
disp('S sparse matrix:')
toc

clear S_cells



%% Q MATRIX

% Approximate number of nonzero elements:

N_nz = round( N_can_add^(7/4) );

% Initialize variable:

Q_cells = zeros(2,N_nz,'uint32');

% Initialize counters:

next = 0;
at_idx = 0;
at0_idx = 0;

% Loop to find desired alters for each type in each preference class
% and check the desire is mutual:

tic
for a = 1:size(Z_names,1)
	z = Z_names(a,:);
    % All pref. classes for this Z:
	H_Z = evalin('caller', ['H_' z]);
	for h = 1:max(H_Z(:,1) )
        % One pref. class:
		H = H_Z(H_Z(:,1) == h, 2);
        can_add_idx = find(can_add(H));
        % Loop over types that can add projects:
        for b = can_add_idx'
			at_idx = at0_idx + b;
            t = H(b);
            % Extract desired alter types:
            temp = Add_cell(t,H);
			s_want = unique( cell2mat( temp(~cellfun('isempty', temp) )' ) );
            if ~isempty(s_want)
                % Positions of desired alter types in list of alloc. par's:
                as_idx = sort( cell2mat( tofa_idx(s_want)' ) );
                as_idx = as_idx(as_idx >= at_idx);
                % Use S to see if they want to add a project w/type t:
                temp = find(S(as_idx,t) );
				len = length(temp);
                % Add them to list of index positions for Q matrix:
				Q_cells(1, next + (1:len) ) = repmat(at_idx,1,len);
				Q_cells(2, next + (1:len) ) = as_idx(temp);
				next = next + len;
            end
        end
        at0_idx = at0_idx + size(H,1);
	end
end
disp('Q matrix cells:')
toc

% Initialize variable:

clear global Q
global Q

% Make sparse matrix from list of index positions:

tic
Q = sparse(double(Q_cells(1,1:next)),double(Q_cells(2,1:next)), true, N_alphas,N_alphas);
Q = Q | Q';
disp('Q sparse matrix:')
toc

clear Q_cells



end
